rRNADepletion_s4ter

Setup

library(tidyverse)
library(magrittr)
library(ngsReports)
library(here)
library(scales)
library(ggpubr)
library(kableExtra)
library(AnnotationHub)
library(ensembldb)
library(edgeR)
library(corrplot)
library(DT)
library(pander)
library(ggrepel)
library(pheatmap)
library(ggdendro)
library(Biostrings)
library(Gviz)
library(Rsamtools)

if (interactive()) setwd(here::here())
theme_set(theme_bw())
cores <- 5

Sequence information

ah_Dr <- AnnotationHub() %>%
  subset(species == "Danio rerio") %>%
  subset(rdataclass == "EnsDb")
ensDb_Dr <- ah_Dr[["AH83189"]]
trEns_Dr <- transcripts(ensDb_Dr) %>%
  mcols() %>% 
  as_tibble()
trLen_Dr <- exonsBy(ensDb_Dr, "tx") %>%
  width() %>%
  vapply(sum, integer(1))
geneGcLen_Dr <- trLen_Dr %>%
  enframe() %>%
  set_colnames(c("tx_id", "length")) %>%
  left_join(trEns_Dr) %>%
  group_by(gene_id) %>% 
  summarise(
    aveLen = mean(length),
    maxLen = max(length), 
    aveGc = sum(length * gc_content) / sum(length),
    longestGc = gc_content[which.max(length)[[1]]]
  ) %>%
  mutate(
    aveGc =  aveGc / 100,
    longestGc = longestGc / 100
  )
trGcLen_Dr <- trLen_Dr %>%
  enframe() %>%
  set_colnames(c("tx_id", "length")) %>%
  left_join(trEns_Dr) %>%
  group_by(tx_id) %>% 
  summarise(
    aveLen = mean(length),
    maxLen = max(length), 
    aveGc = sum(length * gc_content) / sum(length),
    longestGc = gc_content[which.max(length)[[1]]]
  ) %>%
  mutate(
    aveGc =  aveGc / 100,
    longestGc = longestGc / 100
  )
genesGR_Dr <- genes(ensDb_Dr)
mcols(genesGR_Dr) <- mcols(genesGR_Dr)[c("gene_id", "gene_name", 
                                         "gene_biotype", "entrezid")]
txGR_Dr <- transcripts(ensDb_Dr)
mcols(txGR_Dr) <- mcols(txGR_Dr)[c("tx_id", "tx_name", 
                                   "tx_biotype", "tx_id_version", "gene_id")]

An EnsDb object was obtained for Ensembl release 101 using the AnnotationHub package. This provided the GC content and length for every gene and transcript in the release. For zebrafish, this consists of 37241 genes and 65905 transcripts.

Raw data

This is a total RNA-seq dataset generated from a 3-way comparison of WT zebrafish (Danio rerio) with heterozygous mutants (psen2^S4Ter/+) and homozygous mutants (psen2^S4Ter/S4Ter). A previous analysis of this dataset identified the possibility of incomplete ribosomal RNA (rRNA) removal. The following analysis involves an investigation into possible reasons for incomplete rRNA removal and any bias this introduces into the data.

Sample information

files <- list.files(
  path = "/hpcfs/users/a1647910/20200902_Psen2S4Ter/0_rawData/FastQC",
  pattern = "zip",
  full.names = TRUE
)

samples <- tibble(
  sample = str_remove(basename(files), "_fastqc.zip"),
  dataset = NA,
  organism = NA
) %>%
  mutate(
    dataset = ifelse(
      str_detect(sample, "Ps2Ex"), "Psen2S4Ter", dataset
    ),
    organism = ifelse(
      str_detect(sample, "Ps2Ex"), "zebrafish", organism
    )
  )
datasets <- samples$dataset %>% 
  unique()

The following analysis involves 12 paired-end samples across 1 dataset(s): Psen2S4Ter.

Library sizes

rawFqc <- files %>%
  FastqcDataList()

data <- grep("Ps2Ex", fqName(rawFqc))
labels <- rawFqc[data] %>%
  fqName() %>%
  str_remove("_6month_07_07_2016_F3") %>%
  str_remove("\\.fastq\\.gz") %>%
  str_remove("Ps2Ex3M1_")
rawLib <- plotReadTotals(rawFqc[data]) +
  labs(subtitle = "Psen2S4Ter") + 
  scale_x_discrete(labels = labels)

The library sizes of the unprocessed dataset(s) range between 27,979,654 and 37,144,975 reads.

rawLib

GC content

rRNA transcripts are known to have high GC content. Therefore, inspecting the GC content of the raw reads is a logical start point for detecting incomplete rRNA removal. A spike in GC content at > 70% is expected if this is the case.

plotly::ggplotly(
  plotGcContent(
    x = rawFqc[data], 
    plotType = "line",
    gcType = "Transcriptome",
    species = "Drerio"
  ) +
    labs(title = "Psen2S4Ter Dataset (D. rerio)") + 
    theme(legend.position="none")
) 

GC content of reads in the dataset. Clear spikes above 70% GC are observed, which is likely due to incomplete rRNA depletion.

Overrepresented sequences

The top 30 overrepresented sequences were analysed using blastn and were found to be predominantly rRNA sequences.

getModule(rawFqc, "Overrep") %>% 
  group_by(Sequence, Possible_Source) %>% 
  summarise(`Found In` = n(), `Highest Percentage` = max(Percentage)) %>% 
  arrange(desc(`Highest Percentage`), desc(`Found In`)) %>% 
  ungroup() %>% 
  dplyr::slice(1:30) %>%
  mutate(`Highest Percentage` = percent_format(0.01)(`Highest Percentage`/100)) %>%
  kable(
    align = "llrr", 
    caption = paste(
      "Top", nrow(.),"Overrepresented sequences.",
      "The number of samples they were found in is shown,",
      "along with the percentage of the most 'contaminated' sample."
    )
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive")
  )

Top 30 Overrepresented sequences. The number of samples they were found in is shown, along with the percentage of the most ‘contaminated’ sample.

Sequence	Possible_Source	Found In	Highest Percentage
GTGGGTTCAGGTAATTAATTTAAAGCTACTTTCGTGTTTGGGCCTCTAGC	No Hit	12	1.72%
CTGGGGGAGCGGCCGCCCCGCGGCGCCCCCTCTCGTTCCCGTCTCCGGAG	No Hit	10	1.69%
CCGCTGTATTACTCAGGCTGCACTGCAGTGTCTATTCACAGGCGCGATCC	No Hit	12	1.32%
GGCCCGGCGCACGTCCAGAGTCGCCGCCGCACACCGCAGCGCATCCCCCC	No Hit	9	1.31%
CTCCTGAAAAGGTTGTATCCTTTGTTAAAGGGGCTGTACCCTCTTTAACT	No Hit	11	1.11%
GGTTCAGGTAATTAATTTAAAGCTACTTTCGTGTTTGGGCCTCTAGCATC	No Hit	12	1.09%
GGGGTGTACGAAGCTGAACTTTTATTCATCTCCCAGACAACCAGCTATTG	No Hit	12	1.07%
GGCCCGGCGCACGTCCAGAGTCGCCGCCGCGCACCGCAGCGCATCCCCCC	No Hit	10	1.03%
CGAGAGGCTCTAGTTGATATACTACGGCGTAAAGGGTGGTTAAGGAACAA	No Hit	12	1.01%
GGGGGAGCGGCCGCCCCGCGGCGCCCCCTCTCGTTCCCGTCTCCGGAGCG	No Hit	9	0.87%
CCTCCTTCAAGTATTGTTTCATGTTACATTTTCGTATATTCTGGGGTAGA	No Hit	12	0.82%
CCCGCTGTATTACTCAGGCTGCACTGCAGTGTCTATTCACAGGCGCGATC	No Hit	12	0.79%
GTTCAGGTAATTAATTTAAAGCTACTTTCGTGTTTGGGCCTCTAGCATCT	No Hit	12	0.79%
GGGTTCAGGTAATTAATTTAAAGCTACTTTCGTGTTTGGGCCTCTAGCAT	No Hit	12	0.78%
CGGGTCGGGTGGGTGGCCGGCATCACCGCGGACCTCGGGCGCCCTTTTGG	No Hit	12	0.73%
GGGCCTCTAGCATCTAAAAGCGTATAACAGTTAAAGGGCCGTTTGGCTTT	No Hit	11	0.68%
CAGTGGCGTGCGCCTGTAATCCAAGCTACTGGGAGGCTGAGGCTGGCGGA	No Hit	11	0.64%
CGGGTCGGGTGGGTAGCCGGCATCACCGCGGACCTCGGGCGCCCTTTTGG	No Hit	12	0.57%
CTTAGACGACCTGGTAGTCCAAGGCTCCCCCAGGAGCACCATATCGATAC	No Hit	11	0.54%
AGCTGGGGAGATCCGCGAGAAGGGCCCGGCGCACGTCCAGAGTCGCCGCC	No Hit	11	0.53%
GGCCTCTAGCATCTAAAAGCGTATAACAGTTAAAGGGCCGTTTGGCTTTA	No Hit	10	0.53%
CAGCCTATTTAACTTAGGGCCAACCCGTCTCTGTGGCAATAGAGTGGGAA	No Hit	12	0.51%
GGGTGGGTGGCCGGCATCACCGCGGACCTCGGGCGCCCTTTTGGACGTGG	No Hit	10	0.50%
GGGAGCGGCCGCCCCGCGGCGCCCCCTCTCGTTCCCGTCTCCGGAGCGCG	No Hit	9	0.49%
CTGGGAGATGAATAAAAGTTCAGCTTCGTACACCCCAAATTAAAAAATTA	No Hit	10	0.48%
GCCTATTTAACTTAGGGCCAACCCGTCTCTGTGGCAATAGAGTGGGAAGA	No Hit	12	0.47%
GGTCGGGTGGGTGGCCGGCATCACCGCGGACCTCGGGCGCCCTTTTGGAC	No Hit	11	0.45%
CCCCCGAACCCTTCCAAGCCGAACCGGAGCCGGTCGCGGCGCACCGCCGA	No Hit	10	0.45%
GTCGGGTGGGTGGCCGGCATCACCGCGGACCTCGGGCGCCCTTTTGGACG	No Hit	10	0.43%
GCCCACTACGACAACGTGTTTTGTAAATTATGATCTTTATTCTCCTGAAA	No Hit	10	0.43%

Trimmed data

Raw libraries were trimmed using cutadapt v1.14 to remove Illumina adapter sequences. Bases with PHRED score < 30, NextSeq-induced polyG runs and reads shorter than 35bp were also removed.

trimFqc <- list.files(
  path = "/hpcfs/users/a1647910/20200902_Psen2S4Ter/1_trimmedData/FastQC",
  pattern = "zip",
  full.names = TRUE
) %>%
  FastqcDataList()
trimStats <- readTotals(rawFqc) %>%
  dplyr::rename(Raw = Total_Sequences) %>%
  left_join(readTotals(trimFqc), by = "Filename") %>%
  dplyr::rename(Trimmed = Total_Sequences) %>%
  mutate(
    Discarded = 1 - Trimmed/Raw,
    Retained = Trimmed / Raw
  )

After trimming of adapters between 4.16% and 5.08% of reads were discarded.

Aligned data

Trimmed reads were:

1. Aligned to rRNA sequences using the BWA-MEM algorithm to estimate the proportion of reads that were of rRNA origin within each sample. BWA-MEM is recommended for high-quality queries of reads ranging from 70bp to 1Mbp as it is faster and more accurate that alternative algorithms BWA-backtrack and BWA-SW.

2. Aligned to the Danio rerio GRCz11 genome (Ensembl release 101) using STAR v2.7.0d and summarised with featureCounts from the Subread v1.5.2 package. These counts were used for all gene-level analysis.

rRNA proportions

rRnaProp <- read.delim(
  "/hpcfs/users/a1647910/20200902_Psen2S4Ter/3_bwa/log/samples.mapped.all", 
  sep = ":", 
  col.names = c("sample", "proportion"), 
  header = FALSE
) %>% 
  mutate(
    sample = str_remove_all(sample, "_6month_F3|[0-9]*_Ps2Ex3M1_|.mapped"),
    sample = basename(sample),
    proportion = proportion/100,
    dataset = "Psen2S4Ter",
    organism = "zebrafish",
    group = str_extract(sample, "(WT|Heter|Hom)")
  ) %>%
  as_tibble()
rRnaProp$dataset %<>%
  factor(levels = c("Psen2S4Ter"))

rRnaProp %>%
  ggplot(aes(sample, proportion)) +
  geom_bar(stat = "identity", position = "dodge") +
  facet_wrap(~dataset, scales = "free_x") +
  scale_y_continuous(labels = percent) +
  labs(x = "Sample", y = "Percent of Total", fill = "Read pair") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

Percentages of each library that align to rRNA sequences with bwa mem.

rRnaProp %>%
  ggplot(aes(group, proportion, fill = group)) +
  geom_boxplot() +
  scale_y_continuous(labels = percent) +
  labs(x = "Genotype", y = "Percent of total RNA", title = "rRNA proportions of each genotype") +
  scale_fill_discrete(
    name = "Genotype"
  )

Gene GC content and length

dgeList <- read_tsv("/hpcfs/users/a1647910/20200902_Psen2S4Ter/2_alignedData/featureCounts/genes.out") %>%
  set_colnames(basename(colnames(.))) %>%
  set_colnames(str_remove(colnames(.), "Aligned.sortedByCoord.out.bam")) %>%
  set_colnames(str_remove(colnames(.), "_6month_F3")) %>%
  set_colnames(str_remove(colnames(.), "[0-9]*_Ps2Ex3M1_")) %>%
  as.data.frame() %>%
  column_to_rownames("Geneid") %>%
  DGEList() %>%
  calcNormFactors()

dgeList$genes <- genesGR_Dr[rownames(dgeList),]
mcols(dgeList$genes) %<>% 
  as.data.frame() %>% 
  left_join(geneGcLen_Dr)
addInfo <- tibble(
  sample = rRnaProp$sample,
  dataset = "Psen2S4Ter",
  organism = "zebrafish",
  rRNA = rRnaProp$proportion
)
dgeList$samples %<>%
  rownames_to_column("rowname") %>%
  mutate(sample = rowname) %>%
  left_join(addInfo) %>%
  column_to_rownames("rowname")
dgeList$samples$filenames <- read_tsv(
  "/hpcfs/users/a1647910/20200902_Psen2S4Ter/2_alignedData/featureCounts/genes.out"
) %>% 
  dplyr::select(str_subset(colnames(.), "Ps2Ex3M1_")) %>%
  colnames()
dgeList$samples$group <- colnames(dgeList) %>%
  str_extract("(WT|Heter|Hom)") %>%
  factor(levels = c("WT", "Heter", "Hom"))

gcInfo <- function(x) {
  x$counts %>%
    as.data.frame() %>%
    rownames_to_column("gene_id") %>%
    as_tibble() %>%
    pivot_longer(
      cols = colnames(x), 
      names_to = "sample", 
      values_to = "counts"
    ) %>%
    dplyr::filter(
      counts > 0
    ) %>%
    left_join(
      geneGcLen_Dr
    ) %>%
    dplyr::select(
      ends_with("id"), sample, counts, aveGc, maxLen
    ) %>%
    split(f = .$sample) %>%
    lapply(
      function(x){
        DataFrame(
          gc = Rle(x$aveGc, x$counts),
          logLen = Rle(log10(x$maxLen), x$counts)
        )
      }
    ) 
}
gcSummary <- function(x) {
  x %>%
    vapply(function(x){
      c(mean(x$gc), sd(x$gc), mean(x$logLen), sd(x$logLen))
    }, numeric(4)
    ) %>%
    t() %>%
    set_colnames(
      c("mn_gc", "sd_gc", "mn_logLen", "sd_logLen")
    ) %>%
    as.data.frame() %>%
    rownames_to_column("sample") %>%
    as_tibble()
}

rle <- gcInfo(dgeList)

sumGc <- gcSummary(rle)

a <- sumGc %>%
  left_join(dgeList$samples) %>%
  ggplot(aes(rRNA, mn_logLen)) +
  geom_point(aes(colour = group), size = 3) +
  geom_smooth(method = "lm") +
  scale_x_continuous(labels = percent) +
  labs(
    x = "rRNA Proportion of Initial Library",
    y = "Mean log(Length)",
    colour = "Genotype"
  ) 
b <- sumGc %>%
  left_join(dgeList$samples) %>%
  ggplot(aes(rRNA, mn_gc)) +
  geom_point(aes(colour = group), size = 3) +
  geom_smooth(method = "lm") +
  scale_y_continuous(labels = percent) +
  scale_x_continuous(labels = percent) +
  labs(
    x = "rRNA Proportion of Initial Library",
    y = "Mean GC Content",
    colour = "Genotype"
  )
ggarrange(
  a, b, ncol = 2, nrow = 1, 
  common.legend = TRUE, legend = "bottom"
) %>%
  annotate_figure("PsenS4Ter Dataset (D. rerio)")

Comparison of residual bias potentially introduced by incomplete rRNA removal. Regression lines are shown along with standard error bands for each comparison.

PCA

genes2keep <- dgeList %>%
  cpm() %>%
  is_greater_than(1) %>%
  rowSums() %>%
  is_weakly_greater_than(6)
dgeFilt <- dgeList[genes2keep,, keep.lib.sizes = FALSE] %>%
  calcNormFactors()

pca <- cpm(dgeFilt, log = TRUE) %>%
  t() %>%
  prcomp()

pcaCor <- pca$x %>%
  as.data.frame() %>%
  rownames_to_column("sample") %>%
  left_join(sumGc) %>%
  as_tibble() %>% 
  left_join(dgeList$samples) %>%
  dplyr::select(
    PC1, PC2, PC3, 
    Mean_GC = mn_gc, 
    Mean_Length = mn_logLen, 
    rRna_Proportion = rRNA
  ) %>% 
  cor()

a <- pca$x %>%
  as.data.frame() %>%
  rownames_to_column("sample") %>%
  left_join(dgeList$samples) %>%
  as_tibble() %>%
  ggplot(aes(PC1, PC2)) +
  geom_point(aes(colour = group), size = 2) +
  labs(
    x = paste0("PC1 (", percent(summary(pca)$importance["Proportion of Variance","PC1"]),")"),
    y = paste0("PC2 (", percent(summary(pca)$importance["Proportion of Variance","PC2"]),")"),
    colour = "Genotype"
  )
b <- pca$x %>%
  as.data.frame() %>%
  rownames_to_column("sample") %>%
  left_join(dgeList$samples) %>%
  ggplot(aes(PC1, rRNA, label = rRNA)) +
  geom_point(aes(colour = group), size = 2) +
  geom_smooth(method = "lm") +
  geom_text_repel(show.legend = FALSE) +
  scale_y_continuous(labels = percent) +
  labs(
    x = paste0("PC1 (", percent(summary(pca)$importance["Proportion of Variance","PC1"]),")"),
    y = "rRNA Proportion of Initial Library",
    colour = "Genotype"
  )
c <- pca$x %>%
  as.data.frame() %>%
  rownames_to_column("sample") %>%
  left_join(sumGc) %>%
  left_join(dgeList$samples) %>%
  as_tibble() %>%
  ggplot(aes(PC1, mn_logLen)) +
  geom_point(aes(colour = group), size = 2) +
  geom_smooth(method = "lm") +
  labs(
    x = paste0("PC1 (", percent(summary(pca)$importance["Proportion of Variance","PC1"]),")"),
    y = "Mean log(Length)",
    colour = "Genotype"
  )
d <- pca$x %>%
  as.data.frame() %>%
  rownames_to_column("sample") %>%
  left_join(sumGc) %>%
  left_join(dgeList$samples) %>%
  as_tibble() %>%
  ggplot(aes(PC1, mn_gc)) +
  geom_point(aes(colour = group), size = 2) +
  geom_smooth(method = "lm") +
  scale_y_continuous(labels = percent) +
  labs(
    x = paste0("PC1 (", percent(summary(pca)$importance["Proportion of Variance","PC1"]),")"),
    y = "Mean GC",
    colour = "Genotype"
  )
ggarrange(
  a, b, c, d, ncol = 2, nrow = 2,
  common.legend = TRUE, legend = "bottom"
) %>%
  annotate_figure("Psen2S4Ter")

PCA plot showing rRNA proportion, mean GC content and mean log(length) after summarisation to gene-level.

corrplot(
  pcaCor,
  type = "lower", 
  diag = FALSE, 
  addCoef.col = 1, addCoefasPercent = TRUE
)

Correlations between the first three principal components and measured variables: mean GC content, mean log(length) and rRNA proportion.

Differential expression

dgeWTvHet <- dgeFilt[,str_detect(colnames(dgeFilt), "Het|WT")]

design <- model.matrix(~rRNA, data = dgeWTvHet$samples)

voom <- voom(dgeWTvHet, design = design)

fit <- lmFit(voom, design = design)

eBayes <- eBayes(fit)

topTable <- eBayes %>%
  topTable(coef = colnames(design)[2], sort.by = "p", n = Inf) %>%
  set_colnames(str_remove(colnames(.), "ID\\.")) %>%
  mutate(Bonf = p.adjust(P.Value, "bonferroni")) %>%
  mutate(DE = Bonf < 0.05) %>%
  unite(Location, c(seqnames, start, end, width, strand), sep = ":") %>%
  dplyr::select(
    Geneid = gene_id,
    Symbol = gene_name,
    AveExpr,
    logFC,
    P.Value,
    FDR = adj.P.Val,
    Location,
    t,
    DE,
    everything(),
    -B
  ) %>%
  as_tibble()

topTable %>% 
  dplyr::select(Geneid, Symbol, AveExpr, logFC, P.Value, FDR, DE) %>%
  mutate(
    AveExpr = format(round(AveExpr, 2), nsmall = 2),
    logFC = format(round(logFC, 2), nsmall = 2),
    P.Value = sprintf("%.2e", P.Value),
    FDR = sprintf("%.2e", FDR)
  ) %>%
  dplyr::slice(1:200) %>%
  datatable(
    options = list(pageLength = 20), 
    class = "striped hover condensed responsive", 
    filter = "top",
    caption = paste(
      "The top 100 differentially expressed genes.",
      nrow(dplyr::filter(topTable, DE)),
      "of",
      nrow(topTable),
      "genes were classified as DE with an Bonferroni p-value < 0.05."
    )
  )

k-mer analysis

jellyfish v2.3.0 was used to count kmers of trimmed fastq files that had been filtered for rRNA sequences. This was performed for 5 values: \(k = 5, 6, 7, 8, 9, 10\). Lower values of \(k\) lose specificity in comparison to higher values, however as \(k\) increases, the exponential increase of possible kmers causes limitations due to computational processing time.

Setup

k = 5

Counts

k5files <- list.files("/hpcfs/users/a1647910/20200902_Psen2S4Ter/6_jellyfish2pass/k5", pattern = "_dumps.txt", full.names = TRUE)
k5counts <- lapply(k5files, function(x){
  read_delim(x, col_names = c("mer", basename(x)), delim = " ") %>%
    set_colnames(str_remove_all(colnames(.), "_6month_F3|[0-9]*_Ps2Ex3M1_|_dumps\\.txt"))
}) %>%
  purrr::reduce(full_join) %>%
  dplyr::select(mer, contains(c("WT", "Heter")))

k5dge <- k5counts %>%
  as.data.frame() %>%
  column_to_rownames("mer") %>%
  DGEList() %>%
  calcNormFactors()
k5dge$samples %<>%
  rownames_to_column("rowname") %>%
  mutate(sample = rowname) %>%
  left_join(addInfo) %>%
  column_to_rownames("rowname")
k5dge$samples$group <- colnames(k5dge) %>%
  str_extract("(WT|Heter)") %>%
  factor(levels = c("WT", "Heter"))

Properties

k5dist <- k5dge %>%
  cpm(log = TRUE) %>%
  as.data.frame() %>%
  pivot_longer(everything(), names_to = "sample", values_to = "count") %>%
  ggplot(aes(x=count, colour = sample)) +
  geom_density() +
  labs(x = "intensity", title = "Distribution of 5-mers")

k5labels <- k5dge$samples %>% 
  mutate(label = paste0(sample, "\n", percent(rRNA, accuracy = 0.01), " rRNA")) %>% 
  .$label
k5heat <- k5dge %>% 
  cpm(log = TRUE) %>%
  as.data.frame() %>%
  t() %>%
  pheatmap(silent = TRUE, cluster_cols = FALSE,
           show_colnames = FALSE, fontsize = 9,
           fontsize_row = 10, border_color = NA,
           main = "5-mer counts heatmap", labels_row = k5labels)
k5heat$tree_row$labels <- k5labels
k5den <- ggdendrogram(k5heat$tree_row, rotate = TRUE) +
  labs(title = "Hierarchical clustering of 5-mer counts") +
  theme(plot.title = element_text(size = 12))

# Assess cpm values to make sure PCA results are not heavily skewed by highly expressed genes
k5pca <- k5dge %>%
  cpm(log = TRUE) %>%
  t() %>%
  prcomp()

# Quick inspection to check whether first two PCA components capture most of the variability
summary(k5pca)$importance %>% pander(split.tables = Inf)

k5pcaPlot <- k5pca$x %>%
  as.data.frame() %>%
  rownames_to_column("sample") %>%
  as_tibble() %>%
  dplyr::select(sample, PC1, PC2) %>%
  left_join(k5dge$samples) %>%
  ggplot(aes(PC1, PC2, colour = group, label = rRNA)) +
  geom_point(alpha = 0.8, size = 3) +
  geom_text_repel(show.legend = FALSE) +
  labs(
    x = paste0("PC1 (", percent(summary(k5pca)$importance[2, "PC1"]), ")"),
    y = paste0("PC2 (", percent(summary(k5pca)$importance[2, "PC2"]), ")"),
    colour = "Genotype",
    title = "k = 5"
  ) +
  scale_colour_discrete(labels = c("Wildtype", "Mutant"))

k5pcaRrna <- k5pca$x %>%
  as.data.frame() %>%
  rownames_to_column("sample") %>%
  as_tibble() %>%
  dplyr::select(sample, PC1, PC2) %>%
  left_join(k5dge$samples) %>%
  ggplot(aes(PC1, rRNA, label = rRNA)) +
  geom_point(aes(colour = group), alpha = 0.8, size = 3) +
  geom_text_repel(show.legend = FALSE) +
  geom_smooth(method = "lm") +
  labs(
    x = paste0("PC1 (", percent(summary(k5pca)$importance[2, "PC1"]), ")"),
    y = "rRNA proportion",
    colour = "Genotype",
    title = "k = 5"
  ) +
  scale_y_continuous(labels = percent) +
  scale_colour_discrete(labels = c("Wildtype", "Mutant"))

Differential expression

k5design <- model.matrix(~rRNA, data = k5dge$samples)
k5voom <- voom(k5dge, design = k5design)
k5fit <- lmFit(k5voom, design = k5design)
k5eBayes <- eBayes(k5fit)

k5topTable <- k5eBayes %>%
  topTable(coef = colnames(k5design)[2], sort.by = "p", n = Inf) %>%
  set_colnames(str_remove(colnames(.), "ID\\.")) %>%
  rownames_to_column("mer") %>%
  mutate(BY = p.adjust(P.Value, "BY")) %>%
  mutate(DE = BY < 0.05) %>%
  dplyr::select(
    mer,
    AveExpr,
    logFC,
    P.Value,
    FDR = adj.P.Val,
    BY,
    t,
    DE,
    everything(),
    -B
  ) %>%
  as_tibble()

k = 6

Counts

k6files <- list.files("/hpcfs/users/a1647910/20200902_Psen2S4Ter/6_jellyfish2pass/k6", pattern = "_dumps.txt", full.names = TRUE)
k6counts <- lapply(k6files, function(x){
  read_delim(x, col_names = c("mer", basename(x)), delim = " ") %>%
    set_colnames(str_remove_all(colnames(.), "_6month_F3|[0-9]*_Ps2Ex3M1_|_dumps\\.txt"))
}) %>%
  purrr::reduce(full_join) %>%
  dplyr::select(mer, contains(c("WT", "Heter")))

k6dge <- k6counts %>%
  as.data.frame() %>%
  column_to_rownames("mer") %>%
  DGEList() %>%
  calcNormFactors()
k6dge$samples %<>%
  rownames_to_column("rowname") %>%
  mutate(sample = rowname) %>%
  left_join(addInfo) %>%
  column_to_rownames("rowname")
k6dge$samples$group <- colnames(k6dge) %>%
  str_extract("(WT|Heter)") %>%
  factor(levels = c("WT", "Heter"))

Properties

k6dist <- k6dge %>%
  cpm(log = TRUE) %>%
  as.data.frame() %>%
  pivot_longer(everything(), names_to = "sample", values_to = "count") %>%
  ggplot(aes(x=count, colour = sample)) +
  geom_density() +
  labs(x = "intensity", title = "Distribution of 6-mers")

k6labels <- k6dge$samples %>% 
  mutate(label = paste0(sample, "\n", percent(rRNA, accuracy = 0.01), " rRNA")) %>% 
  .$label
k6heat <- k6dge %>% 
  cpm(log = TRUE) %>%
  as.data.frame() %>%
  t() %>%
  pheatmap(silent = TRUE, cluster_cols = FALSE,
           show_colnames = FALSE, fontsize = 9,
           fontsize_row = 10, border_color = NA,
           main = "6-mer counts heatmap", labels_row = k6labels)
k6heat$tree_row$labels <- k6labels
k6den <- ggdendrogram(k6heat$tree_row, rotate = TRUE) +
  labs(title = "Hierarchical clustering of 6-mer counts") +
  theme(plot.title = element_text(size = 12))

# Assess cpm values to make sure PCA results are not heavily skewed by highly expressed genes
k6pca <- k6dge %>%
  cpm(log = TRUE) %>%
  t() %>%
  prcomp()

# Quick inspection to check whether first two PCA components capture most of the variability
summary(k6pca)$importance %>% pander(split.tables = Inf)

# Plot PCA
k6pcaPlot <- k6pca$x %>%
  as.data.frame() %>%
  rownames_to_column("sample") %>%
  as_tibble() %>%
  dplyr::select(sample, PC1, PC2) %>%
  left_join(k6dge$samples) %>%
  ggplot(aes(PC1, PC2, colour = group, label = rRNA)) +
  geom_point(alpha = 0.8, size = 3) +
  geom_text_repel(show.legend = FALSE) +
  labs(
    x = paste0("PC1 (", percent(summary(k6pca)$importance[2, "PC1"]), ")"),
    y = paste0("PC2 (", percent(summary(k6pca)$importance[2, "PC2"]), ")"),
    colour = "Genotype",
    title = "k = 6"
  ) +
  scale_colour_discrete(labels = c("Wildtype", "Mutant"))

k6pcaRrna <- k6pca$x %>%
  as.data.frame() %>%
  rownames_to_column("sample") %>%
  as_tibble() %>%
  dplyr::select(sample, PC1, PC2) %>%
  left_join(k6dge$samples) %>%
  ggplot(aes(PC1, rRNA, label = rRNA)) +
  geom_point(aes(colour = group), alpha = 0.8, size = 3) +
  geom_text_repel(show.legend = FALSE) +
  geom_smooth(method = "lm") +
  labs(
    x = paste0("PC1 (", percent(summary(k6pca)$importance[2, "PC1"]), ")"),
    y = "rRNA proportion",
    colour = "Genotype",
    title = "k = 6"
  ) +
  scale_y_continuous(labels = percent) +
  scale_colour_discrete(labels = c("Wildtype", "Mutant"))

Differential expression

k6design <- model.matrix(~rRNA, data = k6dge$samples)
k6voom <- voom(k6dge, design = k6design)
k6fit <- lmFit(k6voom, design = k6design)
k6eBayes <- eBayes(k6fit)

k6topTable <- k6eBayes %>%
  topTable(coef = colnames(k6design)[2], sort.by = "p", n = Inf) %>%
  set_colnames(str_remove(colnames(.), "ID\\.")) %>%
  rownames_to_column("mer") %>%
  mutate(BY = p.adjust(P.Value, "BY")) %>%
  mutate(DE = BY < 0.05) %>%
  dplyr::select(
    mer,
    AveExpr,
    logFC,
    P.Value,
    FDR = adj.P.Val,
    BY,
    t,
    DE,
    everything(),
    -B
  ) %>%
  as_tibble()

k = 7

Counts

k7files <- list.files("/hpcfs/users/a1647910/20200902_Psen2S4Ter/6_jellyfish2pass/k7", pattern = "_dumps.txt", full.names = TRUE)
k7counts <- lapply(k7files, function(x){
  read_delim(x, col_names = c("mer", basename(x)), delim = " ") %>%
    set_colnames(str_remove_all(colnames(.), "_6month_F3|[0-9]*_Ps2Ex3M1_|_dumps\\.txt"))
}) %>%
  purrr::reduce(full_join) %>%
  dplyr::select(mer, contains(c("WT", "Heter")))

k7dge <- k7counts %>%
  as.data.frame() %>%
  column_to_rownames("mer") %>%
  DGEList() %>%
  calcNormFactors()
k7dge$samples %<>%
  rownames_to_column("rowname") %>%
  mutate(sample = rowname) %>%
  left_join(addInfo) %>%
  column_to_rownames("rowname")
k7dge$samples$group <- colnames(k7dge) %>%
  str_extract("(WT|Heter)") %>%
  factor(levels = c("WT", "Heter"))

Properties

k7dist <- k7dge %>%
  cpm(log = TRUE) %>%
  as.data.frame() %>%
  pivot_longer(everything(), names_to = "sample", values_to = "count") %>%
  ggplot(aes(x=count, colour = sample)) +
  geom_density() +
  labs(x = "intensity", title = "Distribution of 7-mers")

k7labels <- k7dge$samples %>% 
  mutate(label = paste0(sample, "\n", percent(rRNA, accuracy = 0.01), " rRNA")) %>% 
  .$label
k7heat <- k7dge %>% 
  cpm(log = TRUE) %>%
  as.data.frame() %>%
  t() %>%
  pheatmap(silent = TRUE, cluster_cols = FALSE,
           show_colnames = FALSE, fontsize = 9,
           fontsize_row = 10, border_color = NA,
           main = "7-mer counts heatmap", labels_row = k7labels)
k7heat$tree_row$labels <- k7labels
k7den <- ggdendrogram(k7heat$tree_row, rotate = TRUE) +
  labs(title = "Hierarchical clustering of 7-mer counts") +
  theme(plot.title = element_text(size = 12))

# Assess cpm values to make sure PCA results are not heavily skewed by highly expressed genes
k7pca <- k7dge %>%
  cpm(log = TRUE) %>%
  t() %>%
  prcomp()

# Quick inspection to check whether first two PCA components capture most of the variability
summary(k7pca)$importance %>% pander(split.tables = Inf)

# Plot PCA
k7pcaPlot <- k7pca$x %>%
  as.data.frame() %>%
  rownames_to_column("sample") %>%
  as_tibble() %>%
  dplyr::select(sample, PC1, PC2) %>%
  left_join(k7dge$samples) %>%
  ggplot(aes(PC1, PC2, colour = group, label = rRNA)) +
  geom_point(alpha = 0.8, size = 3) +
  geom_text_repel(show.legend = FALSE) +
  labs(
    x = paste0("PC1 (", percent(summary(k7pca)$importance[2, "PC1"]), ")"),
    y = paste0("PC2 (", percent(summary(k7pca)$importance[2, "PC2"]), ")"),
    colour = "Genotype",
    title = "k = 7"
  ) +
  scale_colour_discrete(labels = c("Wildtype", "Mutant"))

k7pcaRrna <- k7pca$x %>%
  as.data.frame() %>%
  rownames_to_column("sample") %>%
  as_tibble() %>%
  dplyr::select(sample, PC1, PC2) %>%
  left_join(k7dge$samples) %>%
  ggplot(aes(PC1, rRNA, label = rRNA)) +
  geom_point(aes(colour = group), alpha = 0.8, size = 3) +
  geom_text_repel(show.legend = FALSE) +
  geom_smooth(method = "lm") +
  labs(
    x = paste0("PC1 (", percent(summary(k7pca)$importance[2, "PC1"]), ")"),
    y = "rRNA proportion",
    colour = "Genotype",
    title = "k = 7"
  ) +
  scale_y_continuous(labels = percent) +
  scale_colour_discrete(labels = c("Wildtype", "Mutant"))

Differential expression

k7design <- model.matrix(~rRNA, data = k7dge$samples)
k7voom <- voom(k7dge, design = k7design)
k7fit <- lmFit(k7voom, design = k7design)
k7eBayes <- eBayes(k7fit)

k7topTable <- k7eBayes %>%
  topTable(coef = colnames(k7design)[2], sort.by = "p", n = Inf) %>%
  set_colnames(str_remove(colnames(.), "ID\\.")) %>%
  rownames_to_column("mer") %>%
  mutate(BY = p.adjust(P.Value, "BY")) %>%
  mutate(DE = BY < 0.05) %>%
  dplyr::select(
    mer,
    AveExpr,
    logFC,
    P.Value,
    FDR = adj.P.Val,
    BY,
    t,
    DE,
    everything(),
    -B
  ) %>%
  as_tibble()

k = 8

Counts

k8files <- list.files("/hpcfs/users/a1647910/20200902_Psen2S4Ter/6_jellyfish2pass/k8", pattern = "_dumps.txt", full.names = TRUE)
k8counts <- lapply(k8files, function(x){
  read_delim(x, col_names = c("mer", basename(x)), delim = " ") %>%
    set_colnames(str_remove_all(colnames(.), "_6month_F3|[0-9]*_Ps2Ex3M1_|_dumps\\.txt"))
}) %>%
  purrr::reduce(full_join) %>%
  dplyr::select(mer, contains(c("WT", "Heter")))

k8dge <- k8counts %>%
  as.data.frame() %>%
  column_to_rownames("mer") %>%
  DGEList() %>%
  calcNormFactors()
k8dge$samples %<>%
  rownames_to_column("rowname") %>%
  mutate(sample = rowname) %>%
  left_join(addInfo) %>%
  column_to_rownames("rowname")
k8dge$samples$group <- colnames(k8dge) %>%
  str_extract("(WT|Heter)") %>%
  factor(levels = c("WT", "Heter"))

Properties

k8dist <- k8dge %>%
  cpm(log = TRUE) %>%
  as.data.frame() %>%
  pivot_longer(everything(), names_to = "sample", values_to = "count") %>%
  ggplot(aes(x=count, colour = sample)) +
  geom_density() +
  labs(x = "intensity", title = "Distribution of 8-mers")

k8labels <- k8dge$samples %>% 
  mutate(label = paste0(sample, "\n", percent(rRNA, accuracy = 0.01), " rRNA")) %>% 
  .$label
k8heat <- k8dge %>% 
  cpm(log = TRUE) %>%
  as.data.frame() %>%
  t() %>%
  pheatmap(silent = TRUE, cluster_cols = FALSE,
           show_colnames = FALSE, fontsize = 9,
           fontsize_row = 10, border_color = NA,
           main = "8-mer counts heatmap", labels_row = k8labels)
k8heat$tree_row$labels <- k8labels
k8den <- ggdendrogram(k8heat$tree_row, rotate = TRUE) +
  labs(title = "Hierarchical clustering of 8-mer counts") +
  theme(plot.title = element_text(size = 12))

# Assess cpm values to make sure PCA results are not heavily skewed by highly expressed genes
k8pca <- k8dge %>%
  cpm(log = TRUE) %>%
  t() %>%
  prcomp()

# Quick inspection to check whether first two PCA components capture most of the variability
summary(k8pca)$importance %>% pander(split.tables = Inf)

# Plot PCA
k8pcaPlot <- k8pca$x %>%
  as.data.frame() %>%
  rownames_to_column("sample") %>%
  as_tibble() %>%
  dplyr::select(sample, PC1, PC2) %>%
  left_join(k8dge$samples) %>%
  ggplot(aes(PC1, PC2, colour = group, label = rRNA)) +
  geom_point(alpha = 0.8, size = 3) +
  geom_text_repel(show.legend = FALSE) +
  labs(
    x = paste0("PC1 (", percent(summary(k8pca)$importance[2, "PC1"]), ")"),
    y = paste0("PC2 (", percent(summary(k8pca)$importance[2, "PC2"]), ")"),
    colour = "Genotype",
    title = "k = 8"
  ) +
  scale_colour_discrete(labels = c("Wildtype", "Mutant"))

k8pcaRrna <- k8pca$x %>%
  as.data.frame() %>%
  rownames_to_column("sample") %>%
  as_tibble() %>%
  dplyr::select(sample, PC1, PC2) %>%
  left_join(k8dge$samples) %>%
  ggplot(aes(PC1, rRNA, label = rRNA)) +
  geom_point(aes(colour = group), alpha = 0.8, size = 3) +
  geom_text_repel(show.legend = FALSE) +
  geom_smooth(method = "lm") +
  labs(
    x = paste0("PC1 (", percent(summary(k8pca)$importance[2, "PC1"]), ")"),
    y = "rRNA proportion",
    colour = "Genotype",
    title = "k = 8"
  ) +
  scale_y_continuous(labels = percent) +
  scale_colour_discrete(labels = c("Wildtype", "Mutant"))

Differential expression

k8design <- model.matrix(~rRNA, data = k8dge$samples)
k8voom <- voom(k8dge, design = k8design)
k8fit <- lmFit(k8voom, design = k8design)
k8eBayes <- eBayes(k8fit)

k8topTable <- k8eBayes %>%
  topTable(coef = colnames(k8design)[2], sort.by = "p", n = Inf) %>%
  set_colnames(str_remove(colnames(.), "ID\\.")) %>%
  rownames_to_column("mer") %>%
  mutate(BY = p.adjust(P.Value, "BY")) %>%
  mutate(DE = BY < 0.05) %>%
  dplyr::select(
    mer,
    AveExpr,
    logFC,
    P.Value,
    FDR = adj.P.Val,
    BY,
    t,
    DE,
    everything(),
    -B
  ) %>%
  as_tibble()

k = 9

Counts

k9files <- list.files("/hpcfs/users/a1647910/20200902_Psen2S4Ter/6_jellyfish2pass/k9", pattern = "_dumps.txt", full.names = TRUE)
k9counts <- lapply(k9files, function(x){
  read_delim(x, col_names = c("mer", basename(x)), delim = " ") %>%
    set_colnames(str_remove_all(colnames(.), "_6month_F3|[0-9]*_Ps2Ex3M1_|_dumps\\.txt"))
}) %>%
  purrr::reduce(full_join) %>%
  dplyr::select(mer, contains(c("WT", "Heter")))

k9dge <- k9counts %>%
  as.data.frame() %>%
  column_to_rownames("mer") %>%
  DGEList() %>%
  calcNormFactors()
k9dge$samples %<>%
  rownames_to_column("rowname") %>%
  mutate(sample = rowname) %>%
  left_join(addInfo) %>%
  column_to_rownames("rowname")
k9dge$samples$group <- colnames(k9dge) %>%
  str_extract("(WT|Heter)") %>%
  factor(levels = c("WT", "Heter"))

Properties

k9dist <- k9dge %>%
  cpm(log = TRUE) %>%
  as.data.frame() %>%
  pivot_longer(everything(), names_to = "sample", values_to = "count") %>%
  ggplot(aes(x=count, colour = sample)) +
  geom_density() +
  labs(x = "intensity", title = "Distribution of 9-mers")

k9labels <- k9dge$samples %>% 
  mutate(label = paste0(sample, "\n", percent(rRNA, accuracy = 0.01), " rRNA")) %>% 
  .$label
k9heat <- k9dge %>% 
  cpm(log = TRUE) %>%
  as.data.frame() %>%
  t() %>%
  pheatmap(silent = TRUE, cluster_cols = FALSE,
           show_colnames = FALSE, fontsize = 9,
           fontsize_row = 10, border_color = NA,
           main = "9-mer counts heatmap", labels_row = k9labels)
k9heat$tree_row$labels <- k9labels
k9den <- ggdendrogram(k9heat$tree_row, rotate = TRUE) +
  labs(title = "Hierarchical clustering of 9-mer counts") +
  theme(plot.title = element_text(size = 12))

# Assess cpm values to make sure PCA results are not heavily skewed by highly expressed genes
k9pca <- k9dge %>%
  cpm(log = TRUE) %>%
  t() %>%
  prcomp()

# Quick inspection to check whether first two PCA components capture most of the variability
summary(k9pca)$importance %>% pander(split.tables = Inf)

# Plot PCA
k9pcaPlot <- k9pca$x %>%
  as.data.frame() %>%
  rownames_to_column("sample") %>%
  as_tibble() %>%
  dplyr::select(sample, PC1, PC2) %>%
  left_join(k9dge$samples) %>%
  ggplot(aes(PC1, PC2, colour = group, label = rRNA)) +
  geom_point(alpha = 0.8, size = 3) +
  geom_text_repel(show.legend = FALSE) +
  labs(
    x = paste0("PC1 (", percent(summary(k9pca)$importance[2, "PC1"]), ")"),
    y = paste0("PC2 (", percent(summary(k9pca)$importance[2, "PC2"]), ")"),
    colour = "Genotype",
    title = "k = 9"
  ) +
  scale_colour_discrete(labels = c("Wildtype", "Mutant"))

k9pcaRrna <- k9pca$x %>%
  as.data.frame() %>%
  rownames_to_column("sample") %>%
  as_tibble() %>%
  dplyr::select(sample, PC1, PC2) %>%
  left_join(k9dge$samples) %>%
  ggplot(aes(PC1, rRNA, label = rRNA)) +
  geom_point(aes(colour = group), alpha = 0.8, size = 3) +
  geom_text_repel(show.legend = FALSE) +
  geom_smooth(method = "lm") +
  labs(
    x = paste0("PC1 (", percent(summary(k9pca)$importance[2, "PC1"]), ")"),
    y = "rRNA proportion",
    colour = "Genotype",
    title = "k = 9"
  ) +
  scale_y_continuous(labels = percent) +
  scale_colour_discrete(labels = c("Wildtype", "Mutant"))

Differential expression

k9design <- model.matrix(~rRNA, data = k9dge$samples)
k9voom <- voom(k9dge, design = k9design)
k9fit <- lmFit(k9voom, design = k9design)
k9eBayes <- eBayes(k9fit)

k9topTable <- k9eBayes %>%
  topTable(coef = colnames(k9design)[2], sort.by = "p", n = Inf) %>%
  set_colnames(str_remove(colnames(.), "ID\\.")) %>%
  rownames_to_column("mer") %>%
  mutate(BY = p.adjust(P.Value, "BY")) %>%
  mutate(DE = BY < 0.05) %>%
  dplyr::select(
    mer,
    AveExpr,
    logFC,
    P.Value,
    FDR = adj.P.Val,
    BY,
    t,
    DE,
    everything(),
    -B
  ) %>%
  as_tibble()

k = 10

Counts

k10files <- list.files("/hpcfs/users/a1647910/20200902_Psen2S4Ter/6_jellyfish2pass/k10", pattern = "_dumps.txt", full.names = TRUE)
k10counts <- lapply(k10files, function(x){
  read_delim(x, col_names = c("mer", basename(x)), delim = " ") %>%
    set_colnames(str_remove_all(colnames(.), "_6month_F3|[0-9]*_Ps2Ex3M1_|_dumps\\.txt"))
}) %>%
  purrr::reduce(full_join) %>%
  dplyr::select(mer, contains(c("WT", "Heter")))

k10dge <- k10counts %>%
  as.data.frame() %>%
  column_to_rownames("mer") %>%
  DGEList() %>%
  calcNormFactors()
k10dge$samples %<>%
  rownames_to_column("rowname") %>%
  mutate(sample = rowname) %>%
  left_join(addInfo) %>%
  column_to_rownames("rowname")
k10dge$samples$group <- colnames(k10dge) %>%
  str_extract("(WT|Heter)") %>%
  factor(levels = c("WT", "Heter"))

Properties

k10dist <- k10dge %>%
  cpm(log = TRUE) %>%
  as.data.frame() %>%
  pivot_longer(everything(), names_to = "sample", values_to = "count") %>%
  ggplot(aes(x=count, colour = sample)) +
  geom_density() +
  labs(x = "intensity", title = "Distribution of 10-mers")

k10labels <- k10dge$samples %>% 
  mutate(label = paste0(sample, "\n", percent(rRNA, accuracy = 0.01), " rRNA")) %>% 
  .$label
k10heat <- k10dge %>% 
  cpm(log = TRUE) %>%
  as.data.frame() %>%
  t() %>%
  pheatmap(silent = TRUE, cluster_cols = FALSE,
           show_colnames = FALSE, fontsize = 9,
           fontsize_row = 10, border_color = NA,
           main = "10-mer counts heatmap", labels_row = k10labels)
k10heat$tree_row$labels <- k10labels
k10den <- ggdendrogram(k10heat$tree_row, rotate = TRUE) +
  labs(title = "Hierarchical clustering of 10-mer counts") +
  theme(plot.title = element_text(size = 12))

# Assess cpm values to make sure PCA results are not heavily skewed by highly expressed genes
k10pca <- k10dge %>%
  cpm(log = TRUE) %>%
  t() %>%
  prcomp()

# Quick inspection to check whether first two PCA components capture most of the variability
summary(k10pca)$importance %>% pander(split.tables = Inf)

# Plot PCA
k10pcaPlot <- k10pca$x %>%
  as.data.frame() %>%
  rownames_to_column("sample") %>%
  as_tibble() %>%
  dplyr::select(sample, PC1, PC2) %>%
  left_join(k10dge$samples) %>%
  ggplot(aes(PC1, PC2, colour = group, label = rRNA)) +
  geom_point(alpha = 0.8, size = 3) +
  geom_text_repel(show.legend = FALSE) +
  labs(
    x = paste0("PC1 (", percent(summary(k10pca)$importance[2, "PC1"]), ")"),
    y = paste0("PC2 (", percent(summary(k10pca)$importance[2, "PC2"]), ")"),
    colour = "Genotype",
    title = "k = 10"
  ) +
  scale_colour_discrete(labels = c("Wildtype", "Mutant"))

k10pcaRrna <- k10pca$x %>%
  as.data.frame() %>%
  rownames_to_column("sample") %>%
  as_tibble() %>%
  dplyr::select(sample, PC1, PC2) %>%
  left_join(k10dge$samples) %>%
  ggplot(aes(PC1, rRNA, label = rRNA)) +
  geom_point(aes(colour = group), alpha = 0.8, size = 3) +
  geom_text_repel(show.legend = FALSE) +
  geom_smooth(method = "lm") +
  labs(
    x = paste0("PC1 (", percent(summary(k10pca)$importance[2, "PC1"]), ")"),
    y = "rRNA proportion",
    colour = "Genotype",
    title = "k = 10"
  ) +
  scale_y_continuous(labels = percent) +
  scale_colour_discrete(labels = c("Wildtype", "Mutant"))

Differential expression

k10design <- model.matrix(~rRNA, data = k10dge$samples)
k10voom <- voom(k10dge, design = k10design)
k10fit <- lmFit(k10voom, design = k10design)
k10eBayes <- eBayes(k10fit)

k10topTable <- k10eBayes %>%
  topTable(coef = colnames(k10design)[2], sort.by = "p", n = Inf) %>%
  set_colnames(str_remove(colnames(.), "ID\\.")) %>%
  rownames_to_column("mer") %>%
  mutate(BY = p.adjust(P.Value, "BY")) %>%
  mutate(DE = BY < 0.05) %>%
  dplyr::select(
    mer,
    AveExpr,
    logFC,
    P.Value,
    FDR = adj.P.Val,
    BY,
    t,
    DE,
    everything(),
    -B
  ) %>%
  as_tibble()

Distributions

ggarrange(
  k5dist, k6dist, k7dist, k8dist, k9dist, k10dist,
  ncol = 2, nrow = 3, common.legend = TRUE, legend = "bottom" 
)

Distributions of kmer counts for each value of \(k\)

PCA

ggarrange(
  k5pcaPlot, k6pcaPlot, k7pcaPlot, k8pcaPlot, k9pcaPlot, k10pcaPlot,
  ncol = 2, nrow = 3, common.legend = TRUE, legend = "bottom"
)

Principal component analysis for all values of \(k\). As \(k\) increases, WT and mutant groups start to separate.

ggarrange(
  k5pcaRrna, k6pcaRrna, k7pcaRrna, k8pcaRrna, k9pcaRrna, k10pcaRrna,
  ncol = 2, nrow = 3, common.legend = TRUE, legend = "bottom"
)

Principal component analysis for all values of \(k\). As \(k\) increases, WT and mutant groups start to separate.

Clustering

ggarrange(
  k5den, k6den, k7den, k8den, k9den, k10den,
  ncol = 2, nrow = 3, common.legend = TRUE, legend = "bottom"
)

Hierarchical clustering of samples based on kmer counts. As \(k\) increases, WT and mutant groups start to cluster.

Differential expression

topRes <- function(x, cap){
  x %>% 
    dplyr::select(mer, AveExpr, logFC, P.Value, FDR, BY, DE) %>%
    mutate(
      AveExpr = format(round(AveExpr, 2), nsmall = 2),
      logFC = format(round(logFC, 2), nsmall = 2),
      P.Value = sprintf("%.2e", P.Value),
      FDR = sprintf("%.2e", FDR),
      BY = sprintf("%.2e", BY)
    ) %>%
    dplyr::slice(1:100) %>%
    datatable(
      options = list(pageLength = 10),
      class = "striped hover condensed responsive",
      filter = "top",
      caption = cap
    )
}

k = 5

topRes(k5topTable,
       cap = paste(
         "The top 100 differentially expressed 5-mers.",
         nrow(dplyr::filter(k5topTable, DE)),
         "of",
         nrow(k5topTable),
         "detected sequences were classified as DE with BY p-value < 0.05."
       )
)

k5topTable %>%
  ggplot(aes(P.Value)) +
  geom_histogram(
    binwidth = 0.02,
    colour = "black", fill = "grey90"
  ) +
  labs(title = "k = 5")

Histogram of p-values. Values follow the expected distribution when there are many differences.

k5topTable %>%
  ggplot(aes(logFC, -log10(P.Value), colour = DE)) +
  geom_point(alpha = 0.5) +
  scale_colour_manual(values = c("black", "grey50", "red")) +
  geom_text_repel(
    data = . %>%
      # dplyr::filter(-log10(P.Value) > 4 | logFC > 4 | logFC < -2),
      dplyr::filter(-log10(P.Value) > 4 | logFC > 3 | logFC < -2.5),
    aes(label = mer, color = "black")
  ) +
  labs(x = "logFC", y = expression(paste(-log[10], "(p)")), title = "k = 5") +
  theme_bw() +
  theme(legend.position = "none")

Volcano plot showing -log10(p-value) against logFC. Kmers classified as DE with BY p-value < 0.05 are highlighted in red.

---

k = 6

topRes(k6topTable,
       cap = paste(
         "The top 100 differentially expressed 6-mers.",
         nrow(dplyr::filter(k6topTable, DE)),
         "of",
         nrow(k6topTable),
         "detected sequences were classified as DE with a BY p-value < 0.05."
       )
)

k6topTable %>%
  ggplot(aes(P.Value)) +
  geom_histogram(
    binwidth = 0.02,
    colour = "black", fill = "grey90"
  ) +
  labs(title = "k = 6")

Histogram of p-values. Values follow the expected distribution when there are some differences.

k6topTable %>%
  ggplot(aes(logFC, -log10(P.Value), colour = DE)) +
  geom_point(alpha = 0.5) +
  scale_colour_manual(values = c("black", "grey50", "red")) +
  geom_text_repel(
    data = . %>%
      dplyr::filter(-log10(P.Value) > 6 | logFC > 6 | logFC < -2.3),
    aes(label = mer, color = "black")
  ) +
  labs(x = "logFC", y = expression(paste(-log[10], "(p)")), title = "k = 6") +
  theme_bw() +
  theme(legend.position = "none")

Volcano plot showing -log10(p-value) against logFC. Kmers classified as DE with BY p-value < 0.05 are highlighted in red.

---

k = 7

topRes(k7topTable,
       cap = paste(
         "The top 100 differentially expressed 7-mers.",
         nrow(dplyr::filter(k7topTable, DE)),
         "of",
         nrow(k7topTable),
         "detected sequences were classified as DE with a BY p-value < 0.05."
       )
)

k7topTable %>%
  ggplot(aes(P.Value)) +
  geom_histogram(
    binwidth = 0.02,
    colour = "black", fill = "grey90"
  ) +
  labs(title = "k = 7")

Histogram of p-values. Values follow the expected distribution when there are many differences.

k7topTable %>%
  ggplot(aes(logFC, -log10(P.Value), colour = DE)) +
  geom_point(alpha = 0.5) +
  scale_colour_manual(values = c("black", "grey50", "red")) +
  geom_text_repel(
    data = . %>%
      dplyr::filter(-log10(P.Value) > 6.4 | logFC > 10 | logFC < -5),
    aes(label = mer, color = "black")
  ) +
  labs(x = "logFC", y = expression(paste(-log[10], "(p)")), title = "k = 7") +
  theme_bw() +
  theme(legend.position = "none")

Volcano plot showing -log10(p-value) against logFC. Kmers classified as DE with BY p-value < 0.05 are highlighted in red.

---

k = 8

topRes(k8topTable,
       cap = paste(
         "The top 100 differentially expressed 8-mers.",
         nrow(dplyr::filter(k8topTable, DE)),
         "of",
         nrow(k8topTable),
         "detected sequences were classified as DE with a BY p-value < 0.05."
       )
)

k8topTable %>%
  ggplot(aes(P.Value)) +
  geom_histogram(
    binwidth = 0.02,
    colour = "black", fill = "grey90"
  ) +
  labs(title = "k = 8")

Histogram of p-values. Values follow the expected distribution when there are many differences.

k8topTable %>%
  ggplot(aes(logFC, -log10(P.Value), colour = DE)) +
  geom_point(alpha = 0.5) +
  scale_colour_manual(values = c("black", "grey50", "red")) +
  geom_text_repel(
    data = . %>%
      dplyr::filter(-log10(P.Value) > 7.2 | logFC > 12.5 | logFC < -8),
    aes(label = mer, color = "black")
  ) +
  labs(x = "logFC", y = expression(paste(-log[10], "(p)")), title = "k = 8") +
  theme_bw() +
  theme(legend.position = "none")

Volcano plot showing -log10(p-value) against logFC. Kmers classified as DE with BY p-value < 0.05 are highlighted in red.

---

k = 9

topRes(k9topTable,
       cap = paste(
         "The top 100 differentially expressed 9-mers.",
         nrow(dplyr::filter(k9topTable, DE)),
         "of",
         nrow(k9topTable),
         "detected sequences were classified as DE with a BY p-value < 0.05."
       )
)

k9topTable %>%
  ggplot(aes(P.Value)) +
  geom_histogram(
    binwidth = 0.02,
    colour = "black", fill = "grey90"
  ) +
  labs(title = "k = 9")

Histogram of p-values. Values follow the expected distribution when there are many differences.

k9topTable %>%
  ggplot(aes(logFC, -log10(P.Value), colour = DE)) +
  geom_point(alpha = 0.5) +
  scale_colour_manual(values = c("black", "grey50", "red")) +
  geom_text_repel(
    data = . %>%
      dplyr::filter(DE & -log10(P.Value) > 7.5 | logFC < -10 | logFC > 15),
    aes(label = mer, color = "black")
  ) +
  labs(x = "logFC", y = expression(paste(-log[10], "(p)")), title = "k = 9") +
  theme_bw() +
  theme(legend.position = "none")

Volcano plot showing -log10(p-value) against logFC. Kmers classified as DE with BY p-value < 0.05 are highlighted in red.

---

k = 10

topRes(k10topTable,
       cap = paste(
         "The top 100 differentially expressed 10-mers.",
         nrow(dplyr::filter(k10topTable, DE)),
         "of",
         nrow(k10topTable),
         "detected sequences were classified as DE with a BY p-value < 0.05."
       )
)

k10topTable %>%
  ggplot(aes(P.Value)) +
  geom_histogram(
    binwidth = 0.02,
    colour = "black", fill = "grey90"
  ) +
  labs(title = "k = 10")

Histogram of p-values. Values follow the expected distribution when there are many differences.

k10topTable %>%
  ggplot(aes(logFC, -log10(P.Value), colour = DE)) +
  geom_point(alpha = 0.5) +
  scale_colour_manual(values = c("black", "grey50", "red")) +
  geom_text_repel(
    data = . %>%
      dplyr::filter(-log10(P.Value) > 8.5 | logFC > 20 | logFC < -13.5),
    aes(label = mer, color = "black")
  ) +
  labs(x = "logFC", y = expression(paste(-log[10], "(p)")), title = "k = 10") +
  theme_bw() +
  theme(legend.position = "none")

Volcano plot showing -log10(p-value) against logFC. Kmers classified as DE with BY p-value < 0.05 are highlighted in red.

GC exploration

## Calculate GC proportion of kmer
gcKmer <- function(x){
  x %>%
    BString() %>%
    letterFrequency(letters = "GC", as.prob = TRUE) %>%
    unname()
}

# ## This chunk takes about 15 mins to run on 6 threads
# ## Objects are saved as .Rds files to save time 
# k5topTable$mer %>% 
#   mclapply(gcKmer, mc.cores = cores) %>%
#   set_names(k5topTable$mer) %>%
#   unlist() %>%
#   enframe(name = "mer", value = "gcProp") %>%
#   saveRDS("R/files/k5gc.Rds")
# k6topTable$mer %>% 
#   mclapply(gcKmer, mc.cores = cores) %>%
#   set_names(k6topTable$mer) %>%
#   unlist() %>%
#   enframe(name = "mer", value = "gcProp") %>%
#   saveRDS("R/files/k6gc.Rds")
# k7topTable$mer %>% 
#   mclapply(gcKmer, mc.cores = cores) %>%
#   set_names(k7topTable$mer) %>%
#   unlist() %>%
#   enframe(name = "mer", value = "gcProp") %>%
#   saveRDS("R/files/k7gc.Rds")
# k8topTable$mer %>% 
#   mclapply(gcKmer, mc.cores = cores) %>%
#   set_names(k8topTable$mer) %>%
#   unlist() %>%
#   enframe(name = "mer", value = "gcProp") %>%
#   saveRDS("R/files/k8gc.Rds")
# k9topTable$mer %>% 
#   mclapply(gcKmer, mc.cores = cores) %>%
#   set_names(k9topTable$mer) %>%
#   unlist() %>%
#   enframe(name = "mer", value = "gcProp") %>%
#   saveRDS("R/files/k9gc.Rds")
# k10topTable$mer %>% 
#   mclapply(gcKmer, mc.cores = cores) %>%
#   set_names(k10topTable$mer) %>%
#   unlist() %>%
#   enframe(name = "mer", value = "gcProp") %>%
#   saveRDS("R/files/k10gc.Rds")

k5gc <- readRDS("R/files/k5gc.Rds")
k6gc <- readRDS("R/files/k6gc.Rds")
k7gc <- readRDS("R/files/k7gc.Rds")
k8gc <- readRDS("R/files/k8gc.Rds")
k9gc <- readRDS("R/files/k9gc.Rds")
k10gc <- readRDS("R/files/k10gc.Rds")

k8topTable %>%
  left_join(k8gc) %>%
  dplyr::arrange(DE) %>%
  ggplot(aes(gcProp, logFC)) +
  geom_point(aes(colour = DE), alpha = 0.5) +
  scale_colour_manual(values = c("black", "red")) +
  labs(title = "k = 8")

k8topTable %>%
  mutate(stat = -sign(logFC)*log10(P.Value)) %>%
  dplyr::arrange(DE) %>%
  left_join(k8gc) %>%
  ggplot(aes(gcProp, stat)) +
  geom_point(aes(colour = DE), alpha = 0.4) +
  scale_x_continuous(labels = percent) +
  scale_colour_manual(values = c("black", "red")) +
  labs(title = "k = 8")

gcPos <- function(x){
  topT <- get(paste0("k", x, "topTable"))
  gc <- get(paste0("k", x, "gc"))
  topT %>%
    dplyr::filter(DE) %>%
    dplyr::filter(logFC > 0) %>%
    left_join(gc) %>%
    ggplot(aes(gcProp)) +
    geom_density() +
    scale_x_continuous(labels = percent) +
    coord_cartesian(xlim = c(0, 1)) +
    labs(title = paste0("k = ", x, ", logFC > 0"))
}
gcNeg <- function(x){
  topT <- get(paste0("k", x, "topTable"))
  gc <- get(paste0("k", x, "gc"))
  topT %>%
    dplyr::filter(DE) %>%
    dplyr::filter(logFC < 0) %>%
    left_join(gc) %>%
    ggplot(aes(gcProp)) +
    geom_density() +
    scale_x_continuous(labels = percent) +
    coord_cartesian(xlim = c(0, 1)) +
    labs(title = paste0("k = ", x, ", logFC < 0"))
}

a <- gcNeg(7)
b <- gcPos(7)
c <- gcNeg(8)
d <- gcPos(8)
e <- gcNeg(9)
f <- gcPos(9)
g <- gcNeg(10)
h <- gcPos(10)
ggarrange(nrow = 4, ncol = 2, a, b, c, d, e, f, g, h) %>%
  annotate_figure("GC content distribution of DE kmers")

DE overlap

DE sequences

dna <- getGenomeTwoBitFile(ensDb_Dr)

deUp <- topTable %>%
  dplyr::filter(DE) %>%
  dplyr::filter(logFC > 0) %>%
  .$Geneid
upExons <- exonsBy(ensDb_Dr, by = "gene", filter = GeneIdFilter(deUp)) %>%
  lapply(function(x){
    GenomicRanges::reduce(x)
  }) %>%
  GRangesList()
upStrings <- lapply(
  upExons, 
  function(x){
    getSeq(dna, x) %>%
      unlist()
  }
) %>%
  DNAStringSet()
upLengths <- vapply(upStrings, length, numeric(1)) %>%
  enframe(name = "Geneid", value = "Length")

deDown <- topTable %>%
  dplyr::filter(DE) %>%
  dplyr::filter(logFC < 0) %>%
  .$Geneid
downExons <- exonsBy(ensDb_Dr, by = "gene", filter = GeneIdFilter(deDown)) %>%
  lapply(function(x){
    GenomicRanges::reduce(x)
  }) %>%
  GRangesList()
downStrings <- lapply(
  downExons, 
  function(x){
    getSeq(dna, x) %>%
      unlist() 
  }
) %>%
  DNAStringSet()
downLengths <- vapply(downStrings, length, numeric(1)) %>%
  enframe(name = "Geneid", value = "Length")

k = 8

k8deUp <- k8topTable %>%
  dplyr::filter(DE) %>%
  dplyr::filter(logFC > 0) %>%
  .$mer
k8deDown <- k8topTable %>%
  dplyr::filter(DE) %>%
  dplyr::filter(logFC < 0) %>%
  .$mer

k8overlapUU <- mclapply(
  upStrings, 
  function(x){
    vapply(
      k8deUp, 
      function(y){
        countPattern(y, x) %>%
          as.logical() %>%
          sum()
      },
      numeric(1)
    ) %>%
      sum()
  },
  mc.cores = cores
) %>%
  unlist() %>%
  enframe(name = "Geneid", value = "Overlap") %>%
  dplyr::arrange(desc(Overlap)) %>%
  left_join(topTable) %>%
  left_join(upLengths) %>%
  dplyr::select(Geneid, Symbol, gene_biotype, AveExpr, logFC, P.Value, Bonf, Length, Overlap) %>% 
  mutate(Norm = (Overlap / Length) * 1000) %>%
  dplyr::arrange(desc(Norm))

k8overlapDD <- mclapply(
  downStrings, 
  function(x){
    vapply(
      k8deDown, 
      function(y){
        countPattern(y, x) %>%
          as.logical() %>%
          sum()
      },
      numeric(1)
    ) %>%
      sum()
  },
  mc.cores = cores
) %>%
  unlist() %>%
  enframe(name = "Geneid", value = "Overlap") %>%
  left_join(topTable) %>%
  left_join(downLengths) %>%
  dplyr::select(Geneid, Symbol, gene_biotype, AveExpr, logFC, P.Value, Bonf, Length, Overlap) %>% 
  mutate(Norm = (Overlap / Length) * 1000) %>%
  dplyr::arrange(desc(Norm))

k8overlapUD <- mclapply(
  upStrings, 
  function(x){
    vapply(
      k8deDown, 
      function(y){
        countPattern(y, x) %>%
          as.logical() %>%
          sum()
      },
      numeric(1)
    ) %>%
      sum()
  },
  mc.cores = cores
) %>%
  unlist() %>%
  enframe(name = "Geneid", value = "Overlap") %>%
  left_join(topTable) %>%
  left_join(upLengths) %>%
  dplyr::select(Geneid, Symbol, gene_biotype, AveExpr, logFC, P.Value, Bonf, Length, Overlap) %>% 
  mutate(Norm = (Overlap / Length) * 1000) %>%
  dplyr::arrange(desc(Norm))

k8overlapDU <- mclapply(
  downStrings, 
  function(x){
    vapply(
      k8deUp, 
      function(y){
        countPattern(y, x) %>%
          as.logical() %>%
          sum()
      },
      numeric(1)
    ) %>%
      sum()
  },
  mc.cores = cores
) %>%
  unlist() %>%
  enframe(name = "Geneid", value = "Overlap") %>%
  left_join(topTable) %>%
  left_join(downLengths) %>%
  dplyr::select(Geneid, Symbol, gene_biotype, AveExpr, logFC, P.Value, Bonf, Length, Overlap) %>% 
  mutate(Norm = (Overlap / Length) * 1000) %>%
  dplyr::arrange(desc(Norm))

k = 9

k9deUp <- k9topTable %>%
  dplyr::filter(DE) %>%
  dplyr::filter(logFC > 0) %>%
  .$mer
k9deDown <- k9topTable %>%
  dplyr::filter(DE) %>%
  dplyr::filter(logFC < 0) %>%
  .$mer

k9overlapUU <- mclapply(
  upStrings, 
  function(x){
    vapply(
      k9deUp, 
      function(y){
        countPattern(y, x) %>%
          as.logical() %>%
          sum()
      },
      numeric(1)
    ) %>%
      sum()
  },
  mc.cores = cores
) %>%
  unlist() %>%
  enframe(name = "Geneid", value = "Overlap") %>%
  dplyr::arrange(desc(Overlap)) %>%
  left_join(topTable) %>%
  left_join(upLengths) %>%
  dplyr::select(Geneid, Symbol, gene_biotype, AveExpr, logFC, P.Value, Bonf, Length, Overlap) %>% 
  mutate(Norm = (Overlap / Length) * 1000) %>%
  dplyr::arrange(desc(Norm))

k9overlapDD <- mclapply(
  downStrings, 
  function(x){
    vapply(
      k9deDown, 
      function(y){
        countPattern(y, x) %>%
          as.logical() %>%
          sum()
      },
      numeric(1)
    ) %>%
      sum()
  },
  mc.cores = cores
) %>%
  unlist() %>%
  enframe(name = "Geneid", value = "Overlap") %>%
  left_join(topTable) %>%
  left_join(downLengths) %>%
  dplyr::select(Geneid, Symbol, gene_biotype, AveExpr, logFC, P.Value, Bonf, Length, Overlap) %>% 
  mutate(Norm = (Overlap / Length) * 1000) %>%
  dplyr::arrange(desc(Norm))

k9overlapUD <- mclapply(
  upStrings, 
  function(x){
    vapply(
      k9deDown, 
      function(y){
        countPattern(y, x) %>%
          as.logical() %>%
          sum()
      },
      numeric(1)
    ) %>%
      sum()
  },
  mc.cores = cores
) %>%
  unlist() %>%
  enframe(name = "Geneid", value = "Overlap") %>%
  left_join(topTable) %>%
  left_join(upLengths) %>%
  dplyr::select(Geneid, Symbol, gene_biotype, AveExpr, logFC, P.Value, Bonf, Length, Overlap) %>% 
  mutate(Norm = (Overlap / Length) * 1000) %>%
  dplyr::arrange(desc(Norm))

k9overlapDU <- mclapply(
  downStrings, 
  function(x){
    vapply(
      k9deUp, 
      function(y){
        countPattern(y, x) %>%
          as.logical() %>%
          sum()
      },
      numeric(1)
    ) %>%
      sum()
  },
  mc.cores = cores
) %>%
  unlist() %>%
  enframe(name = "Geneid", value = "Overlap") %>%
  left_join(topTable) %>%
  left_join(downLengths) %>%
  dplyr::select(Geneid, Symbol, gene_biotype, AveExpr, logFC, P.Value, Bonf, Length, Overlap) %>% 
  mutate(Norm = (Overlap / Length) * 1000) %>%
  dplyr::arrange(desc(Norm))

Coverage patterns

Setup

options(ucscChromosomeNames = FALSE)

exons_Dr <- exonsBy(
  ensDb_Dr,
  by = "tx",
  columns = c("exon_id", "gene_id", "tx_id", "gene_name"),
  filter = GeneIdFilter(c(deDown, deUp))) %>%
  unlist()
colnames(mcols(exons_Dr)) <- str_remove_all(colnames(mcols(exons_Dr)), "_id") %>%
  str_replace_all("tx", "transcript")
mcols(exons_Dr) <- mcols(exons_Dr)[c("gene_name", "gene", "exon", "transcript")]

order <- dgeList$samples %>%
  dplyr::arrange(filenames) %>%
  as_tibble() %>% 
  rownames_to_column("order") %>% 
  dplyr::arrange(-rRNA)

geneCov <- function(geneSym, kVal){
  ## User input
  sym <- geneSym
  kUp <- get(paste0("k", kVal, "deUp"))
  kDown <- get(paste0("k", kVal, "deDown"))
  
  ## Gene region & genome axis tracks
  goi <- genesGR_Dr[mcols(genesGR_Dr)$gene_name == sym,]$gene_id %>%
    .[1] ## Some symbols have multiple gene IDs, let's take the first
  counts <- dgeList$counts[goi,]
  gm <- subset(exons_Dr, gene == goi) 
  # %>% subset(transcript %in% rownames(dgeFiltTx))
  gmRed <- gm %>% 
    GenomicRanges::reduce()
  st <- GenomicRanges::start(gmRed) %>% 
    min()
  en <- GenomicRanges::end(gmRed) %>% 
    max()
  chr <- seqnames(gmRed)@values
  str <- gmRed %>% 
    strand() %>%
    .@values
  gTrack <- GenomeAxisTrack()
  rrTrack <- GeneRegionTrack(
    gmRed,
    name = "Gene model",
    rotation.title = 0,
    col.title = "black",
    fontface.title = 2,
    fontsize = 14,
    fill = "darkorange"
  )
  rTrack <- GeneRegionTrack(
    gm,
    name = "Tx",
    rotation.title = 0,
    col.title = "black",
    fontface.title = 2,
    fontsize = 14
  )
  
  ## GC content track
  win <- slidingWindows(gmRed, width = 1) %>% 
    unlist()
  start <- win %>%
    start()
  end <- win %>%
    end()
  seq <- getSeq(dna, win)
  gc <- seq %>%
    letterFrequency(letters = "GC", OR = "") %>%
    .[,1]
  tbl <- tibble(
    start = start,
    end = end,
    gc = gc
  ) %>% 
    mutate(
      seqnames = chr,
      strand = str,
      cumsum = cumsum(gc),
      cummean = cummean(gc),
      lag = dplyr::lag(cumsum, n = 50, default = NA),
      lead = dplyr::lead(cumsum, n = 50, default = NA),
      gcsum = lead - lag,
      gcprop = gcsum / 101
    ) 
  gcTrack <- tbl %>% 
    dplyr::select(seqnames, start, end, strand, gcprop) %>%
    makeGRangesFromDataFrame(
      keep.extra.columns = TRUE
    ) %>%
    DataTrack(
      type = c("b"),
      ylim = c(0, 1), 
      cex = 0.2,
      cex.axis = 0.3,
      name = "GC %",
      col.title="black",
      fontface.title=2,
      fontsize=14,
      grid = TRUE,
      lwd.grid = 0.5,
      lty.grid = "dashed",
      v = 0
    )
  
  ## Coverage tracks
  files <- list.files(
    "/hpcfs/users/a1647910/20200902_Psen2S4Ter/2_alignedData/bam", 
    pattern = ".bam$", 
    full.names = TRUE
  )
  param <- ScanBamParam(
    what = c("pos", "qwidth"),
    which = GRanges(seqnames = chr, ranges = IRanges(st, en), strand = str),
    flag=scanBamFlag(isUnmappedQuery=FALSE)
  )
  covFunc <- function(x){
    prop <- order$rRNA[order$filenames == x] %>%
      percent(accuracy = 0.1)
    reads <- counts[order$sample[order$filenames == x]]
    pileup(x, scanBamParam = param) %>%
      dplyr::filter(strand == str) %>%
      mutate(start = pos, end = pos, count = count/reads) %>% 
      dplyr::select(seqnames, start, end, strand, count) %>%
      makeGRangesFromDataFrame(keep.extra.columns = TRUE) %>%
      DataTrack(
        type = c("heatmap"),
        cex = 0.2, 
        name = prop,
        rotation.title = 0,
        col.title = "black",
        fontface.title = 2,
        fontsize = 14,
        showAxis = FALSE
      )
  }
  covList <- mclapply(files, covFunc, mc.cores = cores) 
  covOrd <- covList[as.integer(order$order)]
  
  ## k-mer track
  geneMod <- subset(genesGR_Dr, gene_id == goi)
  geneSeq <- getSeq(dna, geneMod)
  ## Up k-mers
  merLocUp <- lapply(kUp, function(x){
    vmatchPattern(x, geneSeq) %>%
      unlist()
  })
  merStartUp <- rapply(merLocUp, function(x){
    start(x)
  })
  merWidthUp <- rapply(merLocUp, function(x){
    width(x)
  })
  upTrack <- DataTrack(
    name = "Up",
    rotation.title = 0,
    col.title = "black",
    fontface.title = 2,
    fontsize = 14,
    showAxis = FALSE
  )
  if (length(merStartUp) > 0) {
    upTrack %<>%
      HighlightTrack(
        trackList = ., 
        start = start(geneMod) + merStartUp,
        width = merWidthUp,
        chromosome = chr,
        inBackground = TRUE,
        col = "white",
        fill = "seagreen",
        showAxis = FALSE
      )
  }
  ## Down k-mers
  merLocDown <- lapply(kDown, function(x){
    vmatchPattern(x, geneSeq) %>%
      unlist()
  })
  merStartDown <- rapply(merLocDown, function(x){
    start(x)
  })
  merWidthDown <- rapply(merLocDown, function(x){
    width(x)
  })
  downTrack <- DataTrack(
    name = "Down",
    rotation.title = 0,
    col.title = "black",
    fontface.title = 2,
    fontsize = 14,
    showAxis = FALSE
  )
  if (length(merStartDown) > 0){
    downTrack %<>%
      HighlightTrack(
        trackList = .,
        name = "test",
        start = start(geneMod) + merStartDown,
        width = merWidthDown,
        chromosome = chr,
        inBackground = TRUE,
        col = "white",
        fill = "red3"
      )
  }
  
  ## Plot the tracks
  plotTracks(
    c(gTrack, rrTrack, rTrack, covOrd, gcTrack, upTrack, downTrack),
    from = st - 100,
    to = en + 100,
    chromosome = chr,
    sizes = c(1,0.5,length(unique(mcols(gm)$transcript))*0.25,rep(0.5, length(covOrd)),1, 0.5, 0.5), 
    main = sym,
    cex.main = 1.5,
    fontface.main = 1
  )
}

Plots

geneCov("BX649540.1", 8)

geneCov("BX649540.1", 9)

geneCov("BX324177.19", 8)

geneCov("rps23", 8)